Acceleration techniques for regularized Newton 
methods applied to electromagnetic inverse 
medium scattering problems 



Thorsten Hohage*and Stefan Langer' 
March 24, 2010 

Abstract 

We study the construction and updating of spectral preconditioners for 
regularized Newton methods and their application to electromagnetic inverse 
medium scattering problems. Moreover, we show how a Lepskii-type stop- 
ping rule can be implemented efficiently for these methods. In numerical 
examples, the proposed method compares favorably with other iterative reg- 
ularization method in terms of work-precision diagrams for exact data. For 
data perturbed by random noise, the Lepskii-type stopping rule performs 
considerably better than the commonly used discrepancy principle. 



1 Introduction 

In this paper we study the efficient numerical solution of an inverse scattering prob- 
lem for time harmonic electromagnetic waves. The forward problem is essentially 
described by the time-harmonic Maxwell equations 

curl curl E(r) - fi;^n(r)^E(r) = 

for the electric field E. Our aim is to reconstruct a local inhomogeneity of the 
refractive index n of a medium, given far field measurements for many incident 
waves. A more detailed discussion of the forward problem is given in ^ 
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After discretization the inverse problem is described by a nonlinear, ill-conditioned 
system of equations F(x) = y with a function F : D(F) C M.'^^ — t- M^, which is 
infinitely smooth on the subset D(F) C M.'^^ where it is defined. Since the system is 
highly ill-conditioned, we have consider the effects of data noise. Here we assume 
an additive noise model for the observe data y°^^: 

y°^^ = F(x) + e (1) 

The noise vector e is assumed to be a vector of random variables with known finite 
covariance matrix and a known bound on the expectation ||]Ee|| < S. 
In this article we contribute to preconditioning techniques for the Levenberg-Marquardt 
algorithm and the iteratively regularized Gauss-Newton method (IRGNM). These 
methods are obtained by applying Tikhonov regularization with some an initial 
guess bfc and a regularization parameter 7^ to the Newton equations A^hk = 
yobs _ F(xfc). Here := F'[xfc] G M^x^ denotes the Jacobian of F at x^. This 
leads to normal equations of the form 

GjG.h, = Gjg, (2) 

with 

G.:^(^,)«<™ and ( ) . 

The choice = corresponds to the Levenberg-Marquardt algorithm and the 
choice bfc = xq — x^ to the IRGNM. As opposed to the Levenberg-Marquardt 
algorithm as used in optimization we simply choose the regularization parameter 7^ 
of the form 

7k = with 7 > 1. (3) 

Convergence and convergence rates of the IRGNM in an infinite dimensional setting 
have been studied first in [21 [71 US]- For further references and results including a 
convergence analysis of Levenberg-Marquardt algorithm we refer to the monographs 

As an alternative, Hanke [I7| suggested to apply the conjugate gradient (CG) 
method the normal equation AjAfch^ = Aj(y°^'^ — F(xfc)) and use the regularizing 
properties of the CG method applied to the normal equation with early stopping. 
This is referred to as Newton-CG method. Regularized Newton methods with inner 
iterative regularization methods have also been studied by Rieder [3T1 [32]. Finally, 
applying a gradient method to the functional x ^ ^||F(x) — y°'^*^||2 leads to the 
nonlinear Landweber iteration x^+i := x^ — juAj(F(xfc) -y"^*^) first studied in [18]. 
For an overview on iterative regularization methods for nonlinear ill-posed problems 
we refer to [22] . 

A continuation method for inverse electromagnetic medium scattering problems 
with multi-frequency data has been studied in [3]. For an overview on level set 
methods for inverse scattering problems we refer to [TTl [T2] . 
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For the inverse electromagnetic scattering problem studied in this paper the eval- 
uation of F and one row of its Jacobian Ak is very expensive and involves the so- 
lution of a three-dimensional forward scattering problem for many incident waves. 
Therefore, a computation of the full Jacobian is not reasonable, and regularization 
method for the inverse problem should only access via matrix-vector multi- 
plications V I— 7- AfcV and g t— t- Ajg. Hence, from the methods discussed above 
only Landweber iteration and Newton-CG can be implemented directly. However, 
the convergence of Landweber iteration is known to be very slow, which is con- 
firmed by our numerical experiments reported in ^ Preconditioning techniques 
for Landweber iteration have been studied in [13], but it is not clear how to ap- 
ply these techniques to inverse electromagnetic medium scattering problems since 
the operator does not act in Hilbert scales. To use the IRGNM and Levenberg- 
Marquardt, we have to solve the system of equations ^ by iterative methods. It 
turns out that standard iterative solvers need many iterations since the systems 
becomes very ill-conditioned as 7^ — )■ 0. 

For the efficient solution of these linear systems we apply the CG-method and exploit 
its close connection to Lanczos' method. The latter method is used to approximately 
compute eigenpairs of G^Gfc to construct a spectral preconditioner for the CG- 
method. Since the eigenvalues \i > ■ ■ ■ > Xm of AjAk decay at an exponential 
rate, it turns out that the approximations determined by Lanczos' method are well 
suited to construct an efficient spectral preconditioner. Spectral preconditioning 
is reviewed in ^ In ^ we describe how the original method proposed in |20] 
can be improved by the construction of updates of the preconditioner during the 
Newton iteration. For a convergence analysis of the IRGNM in combination with 
the discrepancy principle Q discussed below we refer to |2H [26]. It should be 
mentioned that all known convergence results need some condition restricting the 
degree of nonlinearity of F, and unfortunately none of these conditions has been 
verified for the electromagnetic medium scattering problem. 

An essential element of any iterative regularization method for an ill-posed problem 
is a data-driven choice of the stopping index. The most common rule is Morozov's 
discrepancy principle [3D], which consists in stopping the Newton iteration at the 
first index K satisfying 

||F(xa-) - y°^i <tS < ||F(xfc) - y°^i, 0<k<K, t > 1. (4) 

The discrepancy principle is also frequently used for random noise setting 6 = 
A/E]|e|p. However, it is easy to see that this cannot give good results in the limit 
N 00 (see e.g. [5]), and this is confirmed in our numerical experiments. In 
^ we show how a Lepskii-type stopping rule can be implemented efficiently in 
combination with the regularization method studied in Q 

Finally, in ^we report on some numerical experiments to demonstrate the efficiency 
of the methods proposed in this paper. 



3 



2 electromagnetic medium scattering problem 



The propagation of time-harmonic electromagnetic waves in an inhomogeneous, 
non-magnetic, isotropic medium without free charges is described by the time- 
harmonic Maxwell equations 

curl curl E(r) - k^{1 - a(r))E(r) = 0, r G (5a) 

(see [8]). Here E : — t- C'^ describes the space-dependent part of a time-harmonic 
electromagnetic field of the form (E(r)e~*'^*) with angular frequency u > 0. More- 
over, K, := y/EofiQU denotes the wave number, Eq the electric permittivity of vacuum, 
and the magnetic permeability of vacuum. The refractive index of the medium 
given by 

n(r) = a/1 — a(r) 

is assumed to be C^'"-smooth, real and positive in this paper. Moreover, we assume 
that suppa C -Bi = {r G M"^ : |r| < 1}. Now, given a plane incident wave 

E'(r) = E'(r; d, p) = exp(— mr ■ d)p 

with direction d G S"^ and polarization p G such that p • d = 0, the forward 
scattering problem consists in finding a total field E : — )■ satisfying (5a) such 
that the scattered field E*^ := E — E' satisfies the Silver-Miiller radiation condition 

lim (curl E''(r) x r - m|r|E"(r)) = (5b) 

|r|— >oo 

uniformly for all directions r = r/|r| G S*^. The latter condition implies that E*^ has 
the asymptotic behavior 

E^(r; d, p) = ^E-(f ; d, p) + O J , |r| ^ oo 

with a function E°°(-; d, p) : S*^ — > C'^ called the far field pattern of E*^. It satisfies 
E°°(r;d,p) ■ r = 0. 

The inverse problem studied in this paper is to reconstruct a given measurements 
of E°°(r; d, p) for all f , d G 5^ and p G such that d ■ p = 0. 
The forward scattering problem has an equivalent formulation in terms of the elec- 
tromagnetic Lippmann-Schwinger equation 

E(r) + K^ f $(r-f)a(f)E(f)df + grad f $(r - f) ^^^^ ^j'^^ ■ E(f) dr = E^(r) 
Jbi J Br 1 - «(r) 

(6) 

for r G with the scalar fundamental solution $(r) := exp(iK|r|)/(47r|r|). For the 
numerical solution of the forward scattering problems we use a fast solver of (|6|, 
which converges super-linearly for smooth refractive indices (see 
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We typically use between 3 ■ 32^ = 98304 and 3 ■ 64^ = 786 432 degrees of freedom to 
represent E(-; d, p) for each d, p G S"^. The unknown perturbation a of the refractive 
index is represented by a set of coefficients x G M*'^ with 500 < M < 2 000 using 
tensor products of splines in radial direction and spherical harmonics in angular 
direction (see [20] )■ Moreover, we use 25 incident waves with random incident 
directions d^ and random polarizations pj where the directions d^ were drawn from 
the uniform distribution on S^. The exact data are given by complex numbers 
E°°(— d;; dj, pj) ■ pi for j G {1, • • . , 25} and / G {1, . . . , 100} where the dj and Pj 
were generated in the same way as the dj and pj. This yields a real data vector 
y G of size iV = 2 ■ 25 ■ 100 = 5000. 



3 spectral preconditioning 

3.1 CG method and Lanczos' method 

Let us start by recalling the preconditioned conjugate gradient (CG) method and its 
connection to Lanczos' method (see e.g. [TOl [151 [33]). We consider a preconditioned 
equation 

M-^G^Gh = M-^G^g, (7) 

where G G M^^^ is an arbitrary matrix of rank M, and M G is a symmetric 

and positive definite preconditioning matrix. Although the matrix M~^G^G is not 
symmetric in general, the induced linear mapping in is symmetric and positive 
definite with respect to the scalar product {x,y)j^ := {Mx,y) since 

(M-iG^Gx,y>^ = (x,G"^Gy) = (x,M-iG^Gy>^, 
(M-iG^Gx,x)j^ = (x,G^Gx) = IIGxf > for x ^ 0. 

Therefore, the CG-method applied to ([T]) can be coded as follows: 
Algorithm 1 (Preconditioned conjugate gradient method) 

hO = 0; dO = g; r° = G^d"; pi = = M-^r"; I = 0; 

while ||r'|| > e7||h^|| 

1 = 1 + 1; 
q' = Gp'; 

a, = (r'-i,z'-i>/||q'f; 
h' = h'~i + a^p'; 
d' = d^-l -a^q'; 
r' = G^d'; 
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A = <r',z'>/(r'-i,z'-i>; 
p'+i = + /3jp^ 

The stopping criterion ||r'|| > £:7||h'|| ensures a relative accuracy of e/{l — e) of the 
approximate solution if || ((G'^G)"^ || < I/7 (i.e. 7 = 7^ if G = G^, see pH 125]). 
Quantities arising in Algorithm [T] can be used to approximate the largest eigenvalues 
and corresponding eigenvectors of M~^G'''G as follows: Multiplying z-' = p^'^^—fSp^ 
from the left by llz-' Hj^^G and using the definitions and identities 



q 



:J ■- 



q 



q^ 



q^ 



llzJ 



|m y/Oij+l 



|m v^O^ 



yields 



Gz 



;;0 



Gz^ 



j = l,. ..,/-!. 



The identity a^c^ = d-' ^ — d-' multiplied from the left by (Hq-' Haj) ^G^ together 
with 



M 



I q-' II ctj 



||z^||m 

Iq-'lla,- 




yields 



M^G ' q^ 



-z^-i - 




Putting ([s]) and (|9| together we have for all j = 1, — 1 



M-^G^Gz° = 



M-^G^Gz^' = 



~o V^~i 
-z z , 



These formulas can be rewritten as 

M-^G^GZ; = Z;Tz 



0,...,0,^z^ 



where := (z°, 



, z' ^) and 



J_ I ft 

Ql CK2 Oi\ OL1 



\ 



(9) 



(10) 



If we denote by > . . . > > and vi, . . . , the eigenvalues with corresponding 



eigenvectors of the symmetric and positive definite matrix T/, (10) imphes 

ZTm-^G^GZ^v,- = 0,v„ J = 1, . . . , /. 

Hence, in the case that z' vanishes 6i > . . . > 6i are exact eigenvalues of M~^G^G 
with corresponding eigenvectors Z^vi, . . . , Z^v^. In the typical case 7^ the vec- 
tors Z/Vj usually converge rapidly to the eigenvectors corresponding to the outliers 
in the spectrum of M~^G^G (cf. [lOl [15] and the references on the Kaniel-Paige 
theory therein) and Lanczos' method can be interpreted as a particular case of the 
Rayleigh-Ritz method. This connection can be used to interpret the so-called Ritz 
values 61 > . . . > 61 and the Ritz vectors Z;Vi, . . . , Z/V; as approximations to some 
eigenpairs of M~^G^G. 

If WAW^ is an eigendecomposition of the matrix T; with M = I, i.e. W = 
(wi, . . . , W;-) is orthogonal and A = diag(^[^\ . . . , one can prove the equal- 
ity (see m) 

||G^G(Z,-wO - {Zi^i)df\\ = ^"|w40l, ^ = (11) 

where Wj(/) denotes the bottom entry of Wj. This identity can be used to judge 
the accuracy of the Ritz pairs and to decide which of them to use in the spectral 
preconditioner. 

3.2 Spectral preconditioning with Tikhonov regularization 

Assume now that G is of the special form 

A 

with A G M^^^^. Let ui, . . . , um be ortho normal eigenvectors of A'''A, and let 
Ai, . . . , Aa/ be the corresponding eigenvalues. 

Given eigenpairs (Aj,Uj) for j in some non-empty subset J C {1,2, ...,M}, we 
define a spectral preconditioner for G^G = 7I + A^A by 

M:=7l + 5^A,u.(u,)^. 

Its properties are summarized in the following proposition: 
Proposition 2 Assume that rank(A) = M . Then 

a) is symmetric and positive definite, and its inverse is given by 

M-i = ii+vf-^--V.u;. 



G 
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b) M-iQ^Gx = X + Y^jfj Y ^j) = C^GM-^x. 

c) The spectrum of the preconditioned matrix is given by 

a{M-'G^G) = |l + ^ : J ^ j| u {1}, (12) 

and the eigenvalue 1 has multiplicity H^J. 

d) If ^ ^ 1 is an eigenvalue o/M^^G"^G with corresponding eigenvector u, then 
(7(/x — l),u) is an eigenpair of A^A. 

Proof: M is obviously symmetric, and it is positive definite since all its eigen- 
values are > 7 > 0. The formula for the inverse follows from a straightforward 
computation. 

Let U := span{uj : j G J^}. Identifying matrices with their induced linear map- 
pings, we have Ml^^ = G^G|z^ and Mj^^x = jlluJ-, and U and are invariant 
under all the involved linear mappings. Therefore, 

M-^G^GIi, = II,, = G^GM-i|w, 

M-^G^GLx = -G^GU = G^GM-^U . 

7 

Since G^Gx = J^jLiil + ^j) ^j) all x G M.'^'-' , assertion b) follows, c) is 

obtained from b) by inserting the eigenvectors Uj into the formula. 
If (/i, u) is an eigenpair of M~^G^G and /i 7^ 1, it follows that u G U^. Therefore, 
fiu = ^G^Gu = u + ^A^Au. This implies assertion d). □ 

Remark 3 We comment on the assumption rank(A) = M in Proposition For 
the acoustic medium scattering problem injectivity of the continuous Frechet deriva- 
tive -F'[x] has been shown in fWl Prop. 2.2]. For the electromagnetic medium scat- 
tering problem uniqueness proofs for the nonlinear inverse problem (see 0, [7^ ) can 
be modified analogously to show injectivity of F'\x\. It is easy to see that this implies 
injectivity of A if A is a sufficiently accurate discrete approximation of F'[k\ on a 
finite dimensional subspace. 

4 IRGNM with updated spectral preconditioners 

Spectral preconditioning in Newton methods is particularly useful for exponentially 
ill-posed problems such as the electromagnetic inverse medium scattering problem. 
Typically, Lanczos' method approximates outliers in the spectrum well, whereas 
eigenvalues in the bulk of the spectrum are harder to approximate. Frequently 
the more isolated an eigenvalue is, the better the approximation (see |[23j and [TOl 
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Chapter 7]). For exponentially ill-conditioned problems the spectrum of Gjj^Gk,m 
consists of a small number of isolated eigenvalues and a large number of eigenval- 
ues clustering at 7^. If all the large isolated eigenvalues are found and computed 
accurately, spectral preconditioning reduces the condition number significantly. 
Updating the preconditioner may be necessary for the following reasons: 

• If the matrix Gj^Gk,m has multiple isolated eigenvalues, the Lanczos' method 
approximates at most one Ritz pair corresponding to this multiple eigenvalue. 

• During Newton's method the regularization parameter 7^ tends to zero. Hence, 
if we keep the number of known eigenpairs for the construction of the precon- 
ditioner Mfc m fixed, the number of CG-steps will increase rapidly during our 
frozen Newton method (see [23]). 

In the preconditioned Newton iteration we keep the Jacobian = F'(xm) frozen 
for several Newton steps and replace eq. ([T]) by 

M,-LGj^G,,„h = M.-^Gj^g,, (13a) 

where 

G,„:^(^0 a,.d ( ) . (13b) 



Moreover, given some eigenpairs {(A^™\ u^-™'^) : j G J} of A^A^ with orthonormal 
eigenvectors u^"^\ the spectral preconditioner is defined by 

M,,„ := 7.1 + Yl Ar^uJ.'-Vry ■ (14) 

A preconditioned semi-frozen Newton method with updates of the preconditioner 
can be coded as follows: 

Algorithm 4 Input: initial guess xq, data y*^, 5 and/or Gov,, (see ([T 
k := 0; m := 
repeat 



+ Evaluate F(xfc) and define Gk,m,Sk by (13b); 
if Vk + 1 > y/m + 1 + 1 

• m := k; 

• Solve Gj;,Gfc_fchfc = Gj^f^gk by CG-method; 

• Compute via Lanczos' method orthonormal Ritz pairs u^™"*) 
J e J} of G^.G,,,; 
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Select subset J <Z J (see (pLl|) ) and set A^'"'' := 7a:(/^[™'' — 1) for 



j ^ J (see Prop. 

else 

• Define M^,^ by Q; 

if MustUpdate() (see Remark^ below!) 

r-l/Y 

k,m k,m 



J 



^ Solve M-f G^,^G,,^M-:„/^h, = M-,li'Gl^g, by CG-method; 



^ h, := M'JXh,; 

-> Compute Ritz pairs {{fi^f\ uj™^) : j G J2} of M^J^^Gj^Gfc.mM^ 

using Lanczos' method; 
— 7- Select subset J^2 C (^see Remark 6) and set A^'"'' := 7A;(/^j™^ — 

1) for J G J2; 

— )■ Set JT" := JT" U and reorthogonalize {u^-'"'' : j G J"}; 
else 

^ Solve M-,lGl^Gk,mhk = M^j^Gj^^g^ by CG-method; 
end 

end 

+ Xfc+i := Xfc + hk] k := k + 1; 



until StopO (see 
Select stopping index K (see ^ and return xj^; 
We add some remarks on heuristics and implementation details for Algorithm |4} 

1. Usually round-off errors cause loss of orthogonality in the residual vectors 
computed in Algorithm [l| This loss of orthogonality is closely related to the 
convergence of the Ritz vectors (see [I0l[23]). To sustain stability, Algorithm [l] 
was amended by a complete reorthogonalization scheme based on Householder 
transformations (see [T5]). 

2. The necessity of reorthogonalization is also our reason for preconditioning 
with ^ from both sides instead of M^^ from the left when updating the 
preconditioner. In the latter case, reorthogonalization would have to be per- 
formed with respect to the inner product (•, ■)Mfc which is more complicated. 
Note that 



MrL^'x = -^x+ > • I , ^ -^1 (ujx)u 




1 



k,m ^ - pr-^^ z_^\ I — P^\^ i 



\/lk + Aj VTfc^ 
Mlf^x = y^x + ^ (v^TfcTA^ - ^7^) (ujx)uj-. 



10 



3. Spectrally preconditioned linear systems react very sensitively to errors in the 
eigenelements (see [HI El])- Hence, to ensure efficiency of the preconditioner 
it is necessary that the approximations of the Ritz pairs used in the construc- 
tion of the preconditioner be of high accuracy. This is achieved by choosing 
£ = 10~^ in Algorithm [l] when updating or recomputing the preconditioner, 
whereas e = 1/3 is sufficient otherwise. Numerical experience shows that 
computation time invested into improved accuracy of the Ritz pairs pays off 
in the following Newton steps. 

4. MustUpdate(): We update the preconditioner if the last update or recompu- 
tation is at least 4 Newton steps ago and the number of inner iterations in 
the previous Newton step is > 5. 

5. We found it useful not to perform a complete recomputation of the current 
preconditioner if it works well. Therefore, we amend the condition a/A; + 1 > 
\/m + 1 + 1 by the additional requirement that the number of inner iterations 
in the previous step be not too small, say > 8. The condition y/k + 1 > 
■\/m + 1 + 1 is a generalization of the rule to recompute the preconditioner 
whenever /c + 1 is a square number, which was proposed in the original paper 
|20j . Under certain conditions it was shown in j25] to be optimal among all 
rules where is replaced by a function x ^ with /i G (0, 1]. 

6. For updating the preconditioner we only select Ritz values of ^ Gj,! „Gfc_mM 
which are sufficiently well separated from the cluster at 1, say > 1.1. First, 
these eigenvalues are usually computed more accurately by Lanczos' method, 
and second, they are more relevant for preconditioning. 

7. In the initial phase when the updates are large, keeping the Jacobian frozen 
is not efficient. Therefore, we use other methods in this phase, e.g. Newton- 
CG. In some cases globalization strategies will be necessary in this phase, 
although this was not the case in the examples reported below. 

5 Implementation of a Lepskii-type stopping rule 

Lepskii-type stopping rules for regularized Newton methods have been studied in 
lU U] . We refer to the original paper [27] on regression problems and to [2^ I2B] for a 
considerable simplification of the idea and its application to linear inverse problems. 
As opposed to the discrepancy principle, Lepskii-type stopping rules yield order 
optimal rates of convergence for all smoothness classes up to the qualification of 
the underlying linear regularization method (in case of random noise typically only 
up to a logarithmic factor). 

A crucial element of Lepskii's method are estimates of the propagated data noise 
error, and the performance depends essentially on the sharpness of these estimates. 
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Let Rfc := (GjGfc) ^A^. If e G is a deterministic noise vector, an estimate of 
the propagated data noise error is given by 

||Rfce|| < llRfell ||e|| < 7T^||e||, (15) 

and these estimates are sharp if (7fc,e) is an eigenpair of A^Aj . However, if e 
is a random vector with Ee = 0, finite second moments with covariance matrix 



CoVe = (Cov(ej, ej))ij=i..7v, the estimate (15) is usually very pessimistic, and we 
have 



VEllRfceP = Wtrace(R^Cov,Rfc). (16) 



Denoting the right hand side of (15) or (16), respectively, by $(fc), the Lepskii 
stopping rule is defined by 

i^bai := min{A; < i^max : \\^k - x„|| < p$(m), m = k + 1, . . . , Kmax} (17) 

with a parameter p > 4 and a maximal Newton step number -ft'max- We choose 
p = 4.1 in our numerical experiments and -R'max '■= niax{A; G M : $(/c) < R} 
with a reasonable upper bound on the size of propagated data noise in the optimal 
reconstruction. R may be an a-priori known bound ||x — Xo||. However, it is 
advisable to choose a smaller value of R to reduce the number of Newton iterations. 
The final results ^k^^^i do not depend critically on R. 



The main computational challenge in the implementation of the stopping rule (17) 
for random noise is the efficient and accurate computation of $(/c). One possibil- 
ity is to generate L independent copies ei, . . . , e/, of the noise vector and use the 
approximation $(fc) ^ (-^"^ llR-fc^dP)^''^- However, this involves the iterative 
solution of L + l instead of 1 least squares system and leads to a tremendous increase 
of the computational cost. 

With the methods described in the previous sections we can construct approxi- 
mations Hl^^ := y^ -^-r UjwJ of Ri., which allow cheap matrix-vector mul- 
tiplications not involving evaluations of the forward mapping F. This yields the 
approximation 

V 1=1 

Here Wj := AmUj/|| AmUj|| denote the approximated left singular vectors of A^, 
which can be computed directly by an appropriately modified Lanczos method (see 
e.g. [IS])- In the case of white noise, i.e. Cov^ = ct^Iat, the expected value of the 
right hand side is given by the simple expression 



(19) 
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Obviously, equality holds in (19) if (A^)jgj^ is a complete set of eigenvalues of 



A^Am (with multiplicities). Under certain assumptions it has been shown in [6] in 
an infinite dimensional setting that the left hand side can be bounded by a small 
constant times the right hand side uniformly in 7^. if (A^)jgj-^ contains all eigen- 
values > 7fc. Our numerical results in section [6] indicate that this approximation is 
sufficiently accurate. 



exact refractive index 



reconstructed refractive index 




I 



1.3 



1.2 



1.1 



Figure 1: An exact refractive index and its reconstruction by the IRGNM with 
updated preconditioner at iteration 23. The plots show the cube [—1, 1]^, the wave 
number is k = 1. 



6 Numerical results 

As a test example we use the refractive index shown in Figure [ij For further 
information on the forward problem and its numerical solution we refer to |2] and 



Figure |2] illustrates the performance of the IRGNM with updated preconditioners. 
In the update steps for the preconditioner at A; = 16 and 24 the new singular 
values mainly fall into two categories: First, we have singular values which are 
not well separated from the cluster for the regularization parameter 7^ but are 
well separated for 7^.. These singular values are in or near the interval [1/7^, a/ow]- 
The second category are multiple or nearly multiple singular values where only one 
element in the eigenspace is found in the application of Lanczos' method. The 
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stopping rule 


stopping index 


error at stopping index 


optimal 


15.93 ±0.57 


0.0406 ±0.0010 


Lepskii 


12.20 ±0.54 


0.0474 ±0.0016 


discrepancy 


4.87 ±0.43 


0.0744 ± 0.0001 



Table 1: Performance of stopping rules averaged over 15 noise samples for the 
problem in Figure [T] The numbers indicate means and standard deviations. 

use of an update clearly reduces the number of inner CGNE steps in the following 
Newton iterations. 

Moreover, in Fig |3] we compared the speed of convergence of the iterative regu- 
larization methods discussed in the introduction for exact data. Here we measure 
"speed" both in terms of cpu-time and in terms of the number of evaluations of 
F, F'[xm] or F'[xm]'. Landweber iteration is clearly the slowest method although 
some good progress is achieved in the first few steps. The Newton-CG method 
performs very well up to some accuracy after which on it becomes slow, a behavior 
also observed in most other examples. We stopped the Newton-CG iteration at 
an L^-error of ~ 0.28, which was achieved by the updated preconditioned IRGNM 
2.5 times earlier. We also include a comparison with the preconditioned IRGNM 
without updating as suggested in |20]. The updating improves the performance 
particularly at high accuracies. Note that in the first Newton steps where a is still 
small, the iterative solution of the forward problem is faster than in later Newton 
steps. 

Finally, we tested the performance of the Lepskii-type stopping for randomly per- 
turbed data. More precisely, we added independent Gaussian variables to each data 
point. The "relative noise level" ||e||/||y|| was about 2%, but we stress that such a 
point-wise definition of the noise level does not make sense for random noise when 
considering the limit — > oo. We compare the discrepancy principle with r = 2 
to Lepskii's method with p = 4.1. Moreover, we look at the optimal stopping index 
for each noise sample. As expected, the discrepancy principle stops the iteration 
too early. Note in Figure [2] that ||F(xjt) — F(x)|| is at least an order of magnitude 
smaller than at the optimal k 16. (In Figure [2] we used exact data, but the 
behavior is similar for noisy data.) The results in Tab. [l] indicate that Lepskii's 
stopping rule is stable and yields considerably better results than the discrepancy 
principle. 
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\} error of reconstructions 



computed singular values 




Figure 2: Performance of the preconditioned Newton method apphed to the exam- 
ple in Figure [!} The left panels show the continuous L^-error of the reconstructed 
refractive indices and the norm of the residuals ||F(xfc)— y°'^''||2 over the Newton step 
k. The upper right panel shows the computed singular values used for precondi- 
tioning (solid horizontal lines). After step 11 the method changed from Newton-CG 
to IRGNM. At steps m = 11 and 20 the preconditioner was completely recomputed 
for the derivative at a new iterate. This is indicated by solid vertical lines. Updates 
of the preconditioner, indicated by dotted vertical lines, were performed at steps 
16 and 24. The dots on the diagonal line indicate the values of the regularization 
parameter ^7^. Some sufficiently accurate singular values computed during the 
Newton-CG phase are shown for their own interest although they were not used in 
the computation. The right lower panel shows the number of inner CGNE iterations 
over the Newton step k 
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Figure 3: Work-precision diagram comparing the performance of different metliods 
applied to the problem in Figure [T] for exact data. 
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